global SSDIMed "/disk/agedisk4/medicare.work/miller-DUA50377/proj_ssdi"
* Settings
version 16
do "$SSDIMed/scripts/_auxiliary/_project_settings.do"
do "$SSDIMed/scripts/_auxiliary/_project_grstyle.do"

* All required Stata modules are available in the /_packages and /_auxiliary folders
cap adopath - "/home/site/etc/stata/ado.nber"
cap adopath - "/home/site/etc/stata/ado.nber/updates"
adopath ++ "$SSDIMed/scripts/_packages"
adopath ++ "$SSDIMed/scripts/_auxiliary"
grstyle init 





****THIS FILE PREPARES MTO FIGURES WITH CIs for SIX PANELS IN fig: modelrobustness6
***Specs for Figure fig:modelrobustness6 (A.9) that reports marginal medical spending for various entry and spending specs 
***Panel A= 54 & 55
***Panel B= Entry: County. Spending: Years enrolled X County, County X year of observation
***Panel C= Entry: Entry-Month. Spending: Years enrolled X County 
***Panel D= Entry: County. Spending: years enrolled X entry month 
***Panel E= Entry: null. Spending: Years enrolled X county 
***Panel F: Entry: county. Spending: Years enrolled. 


***REPEATS CALCULATIONS OF 04E_MARGINAL_FUNCTIONS_CALCULATIONS.DO

*************************************************************************************************

foreach panel in A B C D E F baseline{
*foreach panel in A{
***describe specs 
if "`panel'"=="A" local entryspec "_spec5657" 
if "`panel'"=="A" local spendingspec "_5657" 
if "`panel'"=="B" local entryspec "_spec02" 
if "`panel'"=="B" local spendingspec "_05"
if "`panel'"=="C" local entryspec "_spec03" 
if "`panel'"=="C" local spendingspec "_02"
if "`panel'"=="D" local entryspec "_spec02" 
if "`panel'"=="D" local spendingspec "_06"
if "`panel'"=="E" local entryspec "_spec01" 
if "`panel'"=="E" local spendingspec "_02"
if "`panel'"=="F" local entryspec "_spec02" 
if "`panel'"=="F" local spendingspec "_01"
if "`panel'"=="baseline" local entryspec "_spec02" 
if "`panel'"=="baseline" local spendingspec "_02"

di "`panel'" 
di "`entryspec'"
di "`spendingspec'"
clear 
use $SSDIMed/results/estimates/x-age52xUR/bootstrap_incidence`entryspec'.dta
compress _all 
rename (intercept age52 UR age52xUR ) (alpha alpha52 alphaUR alpha52XUR)
sum sample
isid sample 
merge 1:1 sample using $SSDIMed/results/estimates/x-age52xUR/bootstrap_spending`spendingspec'.dta, nogen  
d 
if `r(N)'<501 stop  
compress _all 
rename (intercept age52 UR age52xUR ) (beta beta52 betaUR beta52XUR)

***benefit parameters
gen double m=2*beta52/alpha52 
gen double n=beta-(alpha*beta52/alpha52)
gen double mUR=2*(beta52+beta52XUR)/(alpha52+alpha52XUR) 
gen double nUR=beta+betaUR- (alpha+alphaUR)*(beta52+beta52XUR)/(alpha52+alpha52XUR) 
gen double mURminusm=2*(beta52+beta52XUR)/(alpha52+alpha52XUR) - 2*beta52/alpha52
gen double nURminusn=betaUR- (alpha+alphaUR)*(beta52+beta52XUR)/(alpha52+alpha52XUR) + (alpha*beta52/alpha52)
***elements of the figure
gen double x51=alpha 
gen double x51UR=alpha+alphaUR 
gen double x52=alpha+alpha52 
gen double x52UR=alpha+alpha52+alphaUR+alpha52XUR
gen double B_x51=n+m*x51 
gen double B_x51UR=nUR+mUR*x51UR 
gen double B_x52=n+m*x52
gen double B_x52UR=nUR+mUR*x52UR
***difference at x51UR and x52UR
gen double BU_x51=nUR+mUR*x51
gen double BU_x52=nUR+mUR*x52
***set up CI 
gen double B_200=n+m*200
gen double BUR_200=nUR+mUR*200 
gen double B_800=n+m*800
gen double BUR_800=nUR+mUR*800 
if "`panel'"=="A"{
replace B_200=n+m*400
replace BUR_200=nUR+mUR*400 
replace B_800=n+m*1000
replace BUR_800=nUR+mUR*1000 
}
***get pval for comparison of B_x`n' and BU_x`n'


foreach n in 51 52{
sum B_x`n' in 1 
local B_x`n'=`r(mean)'
gen shareBSBU_x`n'belowB_x`n'=(BU_x`n'<`B_x`n'')
replace shareBSBU_x`n'belowB_x`n'=. in 1 
sum BU_x`n' in 1 
local BU_x`n'=`r(mean)'
gen shareBSB_x`n'aboveBU_x`n'=(B_x`n'>`BU_x`n'')
replace shareBSB_x`n'aboveBU_x`n'=. in 1
}
*/

*****COST PARAMETERS 
foreach DeltaC in -1000 -5000{
local posDeltaC=-1*`DeltaC'
gen double mC_51`posDeltaC'=(-`DeltaC' -m*alpha -n +mUR *(alpha+alphaUR) +nUR)/alphaUR
gen double nC_51`posDeltaC'=m*alpha +n- mC_51`posDeltaC'*alpha
gen double alphaCF_51_C`posDeltaC' =(`DeltaC' +nC_51`posDeltaC'-n)/(m-mC_51`posDeltaC')
gen double mC_52`posDeltaC'=(-`DeltaC'  +mUR*(alpha+alphaUR+alpha52+alpha52XUR) -m*(alpha+alpha52) +nUR - n)/(alphaUR+alpha52XUR)
gen double nC_52`posDeltaC'=m*(alpha+alpha52)+n - mC_52`posDeltaC'*(alpha+alpha52)
gen double alphaCF_52_C`posDeltaC' =(`DeltaC' +nC_52`posDeltaC'-n)/(m-mC_52`posDeltaC')
gen double percentCF_health`posDeltaC'=(x51UR -alphaCF_51_C`posDeltaC' + x52UR-alphaCF_52_C`posDeltaC')/(alphaUR+alphaUR+alpha52XUR)
}



*5) collapse and print table 
cap drop ctrl 
save "$SSDIMed/results/temp/robustnessPanel`panel'.dta", replace 
keep in 1 
xpose, clear varname 
save "$SSDIMed/results/temp/tablecolumn`column'_value_Panel`panel'.dta", replace 
use "$SSDIMed/results/temp/robustnessPanel`panel'.dta", clear
drop in 1 
gcollapse (mean) shareBS*, fast 
xpose, clear varname 
save "$SSDIMed/results/temp/tablecolumn`column'_pval_Panel`panel'.dta", replace
use "$SSDIMed/results/temp/robustnessPanel`panel'.dta", clear
drop in 1 
gcollapse (sd) *, fast 
xpose, clear varname 
save "$SSDIMed/results/temp/tablecolumn`column'_sd_Panel`panel'.dta", replace 
use "$SSDIMed/results/temp/robustnessPanel`panel'.dta", clear
drop in 1 
gcollapse (p5) *, fast 
xpose, clear varname 
save "$SSDIMed/results/temp/tablecolumn`column'_lb_Panel`panel'.dta", replace 
use "$SSDIMed/results/temp/robustnessPanel`panel'.dta", clear
drop in 1 
gcollapse (p95) *, fast 
xpose, clear varname 
save "$SSDIMed/results/temp/tablecolumn`column'_ub_Panel`panel'.dta", replace
clear 
gen stat=""
foreach var in value sd lb ub pval{
append using "$SSDIMed/results/temp/tablecolumn`column'_`var'_Panel`panel'.dta"
replace stat="`var'" if stat==""
}
reshape wide v1, i(_varname) j(stat) string 
rename v1* * 
rename _varname var 
drop if var=="sample"
save "$SSDIMed/results/estimates/robustness/tablecolumn`column'_Panel`panel'.dta", replace 

  drop sd 
  drop if (regexm(var, "500") & !regexm(var, "5000")) | regexm(var, "50000")
  keep if regexm(var, "x51") | regexm(var, "x51UR") | regexm(var, "x52") | regexm(var, "x52UR") | regexm(var, "00")
  drop if regexm(var, "C_") | regexm(var, "diff_")
  *replace var=subinstr(var, "`posDeltaC'", "", 1) 
  gen x=.
  gen xlb=.
  gen xub=.
  foreach var in x51 x51UR x52 x52UR {
    replace x=value if var=="`var'"
    replace xlb=lb if var=="`var'"
    replace xub=ub if var=="`var'"
  }
  
  replace x=200 if var=="B_200" | var=="BUR_200"
  replace x=800 if var=="B_800" | var=="BUR_800"
  if "`panel'"=="A"{
  replace x=400 if var=="B_200" | var=="BUR_200"
  replace x=1000 if var=="B_800" | var=="BUR_800" 
  }
  gen B=. 
  gen Blb=.
  gen Bub=.
  foreach var in B_x51 B_x51UR B_x52 B_x52UR B_200 BUR_200 B_800 BUR_800 {
    replace B=value if var=="`var'"
    replace Blb=lb if var=="`var'"
    replace Bub=ub if var=="`var'"
  }

  gen UR=regexm(var, "UR")
  replace var=subinstr(var, "B_", "", 1)
  replace var=subinstr(var, "BUR_", "", 1)
  keep if inlist(var, "x51", "x51UR", "x52", "x52UR", "200", "800")
  foreach var in B Blb Bub {
    replace `var'=`var'/1000
  }

  ***puts everything in one line 
  gcollapse x xlb xub B Blb Bub  , by(var UR) 

  foreach var in x51 x51UR x52 x52UR {
    sum x if var=="`var'" 
    local `var': di %3.0f = `r(mean)'
  }
***below only works when panel==A
 replace var="400" if x==400 
 replace var="1000" if x==1000

**************************************************************************
***DAVID'S FIGURE CODE   
  * Override default graph sizes
  grstyle set size 8pt: text_option
  grstyle set size 9pt: axis_title
  grstyle set size 8pt: tick_label
  grstyle set size 2pt: tickgap
  
  * Arrows and lines
  local p 5
  grstyle set color gs8: p`p'
  grstyle set lpattern solid: p`p'
  grstyle set linewidth vvthin: p`p'
  
  * How many major ticks?
  * How much should range extend beyond ticks (as fraction of tick step size)
  local y_ticks 8
  local y_margin_lo 0.4
  local y_margin_hi 0.4
  
  * y1 tick lo, step size, format
  local y1_tick_lo 11.5
  local y1_step 0.5
  local y1_fmt "%9.1f"
  local y1_sign "$"
  local y1_unit " thousand"
  
  * Construct y labels and scales from parameters above
  forvalues y = 1/1 {
    local y`y'_tick_hi  = `y`y'_tick_lo' + `y`y'_step' * (`y_ticks' - 1)
    local y`y'_range_lo = `y`y'_tick_lo' - `y`y'_step' * `y_margin_lo'
    local y`y'_range_hi = `y`y'_tick_hi' + `y`y'_step' * `y_margin_hi'
    di "y`y'_tick_hi:  `y`y'_tick_hi'"
    di "y`y'_range_lo: `y`y'_range_lo'"
    di "y`y'_range_hi: `y`y'_range_hi'"
    
    local y`y'_label ylabel(`y`y'_tick_lo'(`y`y'_step')`y`y'_tick_hi', axis(`y') format("`y`y'_fmt'") noticks nolabels labgap(1pt)) 
    local y`y'_scale yscale(`yalt' range(`y`y'_range_lo' `y`y'_range_hi') axis(`y') noline) 
    di `"y`y'_label:  `y`y'_label'"'
    di `"y`y'_scale:  `y`y'_scale'"'
  }
  
  * x-axis settings
  local xlo 120
  local xhi 830
  if "`panel'"=="A" local xlo 320
  local xscale xscale(range(`xlo' `xhi'))
  local x_axis xlabel(200(200)800) `xscale' xtitle("Entrants per million", margin(t=-.5))
   if "`panel'"=="A"  local x_axis xlabel(400(200)1000) xscale(range(320 1030)) xtitle("Entrants per million", margin(t=-.5))
   
  * style and color for y-axis title/labels
  local y1_title_style it
  local y1_title_color gs8
  local y1_label_style it
  local y1_label_color gs8
  local y1_label_size 8pt
  
  local text_style sf
  local text_color black
  
  * y-axis labels placed above grid rule (build manually)
  local y1_x `xlo'
  local y1_place ne
  local y1_just left
  local y1_width 25
  
  local y 1
  local y`y'_la
  foreach tick of numlist  `y`y'_tick_lo'(`y`y'_step')`y`y'_tick_hi' {
    if `tick' == `y`y'_tick_hi' {
      * white textbox to go behind label
      local y`y'_la `y`y'_la' text(`=`tick' + `y`y'_step'/100' `y`y'_x' " ", yaxis(`y') place(`y`y'_place') just(`y`y'_just') size(`y`y'_label_size') box lwidth(none) color(`y`y'_label_color') width(`y`y'_width') height(4) bcolor(white%70))
      *label
      local y`y'_la `y`y'_la' text(`tick' `y`y'_x' "{`y`y'_label_style':`y`y'_sign'`=string(`tick', "`y`y'_fmt'")'`y`y'_unit'}", yaxis(`y') place(`y`y'_place') just(`y`y'_just') size(`y`y'_label_size') box lwidth(none) color(`y`y'_label_color') bcolor(none))
    }
    else {
      local y`y'_la `y`y'_la' text(`tick' `y`y'_x' "{`y`y'_label_style':`y`y'_sign'`=string(`tick', "`y`y'_fmt'")'}", yaxis(`y') place(`y`y'_place') just(`y`y'_just') size(`y`y'_label_size') box lwidth(none) color(`y`y'_label_color') bcolor(none))
    }
  }
  
  * y-axis titles, horizontal orientation (build manually)
  local y 1
  local y1title text(`=`y`y'_tick_lo' + `y`y'_step'*(`y_ticks' - 1 + `y_margin_hi')' `xlo' "{`y`y'_title_style':Medical spending}"
  local y1title `y`y'title', yaxis(`y') place(`y`y'_place') just(`y`y'_just') box lwidth(none) color(`y`y'_title_color') bcolor(white) margin(b=+.5))
  local y1title `y`y'title' ytitle(" ", axis(`y') margin(l=-8))
  
  * y-axis settings
  local y_axis1 `y1_label' `y1_la' `y1_scale' `y1title'
  
  * Plot lines
  local y1 B
  local xvar x
  local line_y1_lfit "(lfit B x if UR==0, col(orange) lw(thick) range(200 800))"
  local line_y1_scat "(scatter B x if UR==0 & inrange(x, 201, 799) ,msym(S) mcolor(black) msize(*.8))"
  if "`panel'"=="A" local line_y1_scat "(scatter B x if UR==0 & inrange(x, 401, 999) ,msym(S) mcolor(black) msize(*.8))"
  local line_y2_lfit "(lfit B x if UR==1, lpat(dash) lcol(orange) lw(thick) range(200 800))"
  local line_y2_scat "(scatter B x if UR==1 & inrange(x, 201, 799),msym(Oh) mcolor(black) msize(*1) mlwidth(medthick)) "
    if "`panel'"=="A" local line_y2_scat "(scatter B x if UR==1 & inrange(x, 401, 999),msym(Oh) mcolor(black) msize(*1) mlwidth(medthick)) "
  local ci_area_0 (rarea Blb Bub `xvar' if UR==0, sort col(gs8%30) lw(none))
  local ci_area_1 (rarea Blb Bub `xvar' if UR==1, sort col(gs8%30) lw(none))
  if "`panel'"=="A"   local line_y2_lfit "(lfit B x if UR==1, lpat(dash) lcol(orange) lw(thick) range(400 1000))"
  if "`panel'"=="A"   local line_y1_lfit "(lfit B x if UR==0, col(orange) lw(thick) range(400 1000))"
  
  * Label line_y1
  local text_1 text(12.51 500 "{`text_style':Mean unemployment}"
  local text_1 `text_1', yaxis(1) place(ne) just(left) box lwidth(none) width(45) color(`text_color') bcolor(white%90))
  local xref = `xvar'[4]
  sum `y1' in 4 
  local pci_1 (pci `=r(mean)' `=`xref'' `=r(mean) + .7' `=`xref' + 4.5', pstyle(p5))
  
  * Label line_y2
  local text_2 text(13 550 "{`text_style':+1 p.p. unemployment}"
  local text_2 `text_2', yaxis(1) place(ne) just(left) box lwidth(none) width(45) color(`text_color') bcolor(white%90))
  local xref = `xvar'[4]
  sum `y1' in 4 
  local pci_2 (pci `=r(mean)' `=`xref'' `=r(mean) + 0.07' `=`xref' + 0.45', pstyle(p5))
  
  grstyle set graphsize 3in 3in
  local legend legend(off)
  local region graphregion(margin(t=5 b=-1))
  local ws "placement(w) size(*1)"
  local es "placement(e) size(*1)"
  
  local labels          text(`=B[5]-0.5' `=x[5]-10' "{`text_style':age 49, mean}" "{`text_style':unemployment}", placement(s)) 
  local pci_labels      (pci `=B[5]' `=x[5]' `=B[5]-0.5' `=x[5]-10', pstyle(p5))

  local labels `labels' text(`=B[6]+0.5' `=x[6]+10' "{`text_style':age 49, +1 p.p.}" "{`text_style':unemployment}", placement(n))
  local pci_labels `pci_labels' (pci `=B[6]' `=x[6]' `=B[6]+0.5' `=x[6]+10', pstyle(p5))

  local labels `labels' text(`=B[7]-0.5' `=x[7]-10' "{`text_style':age 50, mean}" "{`text_style':unemployment}", placement(s)) 
  local pci_labels `pci_labels' (pci `=B[7]' `=x[7]' `=B[7]-0.5' `=x[7]-10', pstyle(p5))


  local labels `labels' text(`=B[8]+0.5' `=x[8]+10' "{`text_style':age 50, +1 p.p.}" "{`text_style':unemployment}", placement(n))
  local pci_labels `pci_labels' (pci `=B[8]' `=x[8]' `=B[8]+0.5' `=x[8]+10', pstyle(p5))
if "`panel'"=="A"{
  local labels          text(`=B[5]-0.5' `=x[5]-10' "{`text_style':age 54, mean}" "{`text_style':unemployment}", placement(s)) 
  local pci_labels      (pci `=B[5]' `=x[5]' `=B[5]-0.5' `=x[5]-10', pstyle(p5))

  local labels `labels' text(`=B[6]+0.5' `=x[6]+10' "{`text_style':age 54, +1 p.p.}" "{`text_style':unemployment}", placement(n))
  local pci_labels `pci_labels' (pci `=B[6]' `=x[6]' `=B[6]+0.5' `=x[6]+10', pstyle(p5))

  local labels `labels' text(`=B[7]-0.5' `=x[7]-10' "{`text_style':age 55, mean}" "{`text_style':unemployment}", placement(s)) 
  local pci_labels `pci_labels' (pci `=B[7]' `=x[7]' `=B[7]-0.5' `=x[7]-10', pstyle(p5))


  local labels `labels' text(`=B[8]+0.5' `=x[8]+10' "{`text_style':age 55, +1 p.p.}" "{`text_style':unemployment}", placement(n))
  local pci_labels `pci_labels' (pci `=B[8]' `=x[8]' `=B[8]+0.5' `=x[8]+10', pstyle(p5))
 }
  twoway `ci_area_0' `ci_area_1' `line_y1_lfit' `line_y2_lfit' `line_y1_scat' `line_y2_scat' `pci_labels', `x_axis' `y_axis1' `legend' `region' `labels' graphregion(color(white))  bgcolor(white) 
  
  * Export graph
  cap mkdir "$SSDIMed/results/figures"

  graph export "$SSDIMed/results/figures/robustness_Panel`panel'.pdf", as(pdf) fontface("Source Sans Pro") fontdir("$SSDIMed/data/raw/fonts") replace
}



